clear all; close all; clc;

H = 19100e3;
Rz = 6400e3;
Torb = 12*60*60;
Tmod = Torb/2;
wobr = 2*pi/Torb;
c_light = 3e8;
f_0 =  1602e6;


dt=5*60;
t=0:dt:Tmod;


rr = [0; Rz];
D = nan(1, length(t));
y_rv = nan(1, length(t));

for i = 1:length(t)
    rsv = (H + Rz) * ...
        [cos(wobr*t(i));
        sin(wobr*t(i))]
    %plot([0 rsv(1))], 0 rsv(2)])
    %xlim([-H-RZ; H+Rz])
    %ylim([-Rz; H+Rz])
    %axis equal
    %pause (0,1);
    rv = rsv -rr;
    D(i) = norm(rv);
end
ind = find(y_rv>0);
V = diff(D)/dt;

df = -V/c_light*f_0;
figure;
plot(t, D);
plot(t(2:end), df, t(ind), df(ind-1));
xlabel('t, sec')
ylabel('\Delta f, Hz')

figure;
hist(df)
xlabel('Delta f, Hz')